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On Multi-dimensional Unstructured 
Mesh Adaption 


William A. Wood* and William L. Kleb* 

NASA Langley Research Center, Hampton, VA 23681 

Anisotropic unstructured mesh adaption is developed for a truly multi-dimensional up- 
wind fluctuation splitting scheme, as applied to scalar advection-diffusion. The adaption 
is performed locally using edge swapping, point insertion/deletion, and nodal displace- 
ments. Comparisons are made versus the current state of the art for aggressive anisotropic 
unstructured adaption, which is based on a posteriori error estimates. Demonstration of 
both schemes to model problems, with features representative of compressible gas dy- 
namics, show the present method to be superior to the a posteriori adaption for linear 
advection. The performance of the two methods is more similar when applied to non- 
linear advection, with a difference in the treatment of shocks. The a posteriori adaption 
can excessively cluster points to a shock, while the present multi-dimensional scheme 
tends to merely align with a shock, using fewer nodes. As a consequence of this align- 
ment tendency, an implementation of eigenvalue limiting for the suppression of expansion 
shocks is developed for the multi-dimensional distribution scheme. The differences in the 
treatment of shocks by the adaption schemes, along with the inherently low levels of 
artificial dissipation in the fluctuation splitting solver, suggest the present method is a 
strong candidate for applications to compressible gas dynamics. 


Nomenclature 

A Flux Jacobian 

E r Error estimate 

F Flux function 

l Edge length 

n Outward unit normal 

r Position vector 

S Area of a triangle 

t Time 

U Dependent variable 

a, P Curvilinear advection speeds 

e Eigenvalue limiting parameter 

A Advection velocity 

y Diffusion coefficient 

<j> Fluctuation 

cffi , (jA Fluctuation components 

</>* £ , q N Limited fluctuations 

</> /£ , (/)'’’ Artificial dissipation terms 

<j) v Diffusive fluctuation 

ip Limiter function 

T Minimization functional 

S Weighting matrix 

Cl General element area 

V Gradient operator 

A" Forward difference in iterate n 

Subscripting independent variables represents differ- 
entiation. Tildes indicate averaged states with Prop- 
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erty U. Over-bar indicates cell-averaged value. 

Introduction 

H IGH-FIDELITY aerothermodynamic analyses 
for hypersonic vehicles, such as the X-34, can 
require up to 390 CPU hours on a Cray C-90 for 
a single computational fluid dynamics (CFD) solu- 
tion, even when using a modern, proven, highly-tuned 
solver. 0 The associated solution-adapted grids include 
up to 9 million nodes. It is desired to reduce solu- 
tion times by more than an order of magnitude so 
as to increase the relevancy of CFD to the national 
X-plane programs. Reducing the required mesh sizes, 
while maintaining or improving accuracy, has been tar- 
geted as a high-payoff endeavor for minimizing the 
cost of computational aerothermodynamics. Funda- 
mental algorithmic advances for both the flow solver 
and the mesh-adaption strategy will contribute to re- 
duced costs. 

Solution-adaptive remeshing techniques have been 
utilized with some success for hypersonic flows on 
structured domains. A leader in this field is Gnoffo, 0 
who utilizes a spring analogy energy minimization to 
align the bow shock and cluster to the boundary layer. 
This approach works very well for entry forebodies, for 
which it was developed, but is more difficult to apply 
to complex vehicle shapes. The method also is unre- 
sponsive to embedded shocks or other shock-layer flow 
features. 

Harvey®0’S has developed a mesh adaption tech- 
nique that is sensitive to shock-layer features to obtain 
parabolized Navier-Stokes solutions over simple con- 
figurations, e.g. cones, using a spring analogy based 
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Fig. 1 Pictorial of ‘ideal’ mesh for shock discon- 
tinuity (left) and mesh resulting from curvature- 
based clustering (right). 

on Gnoffo’s work. Unfortunately, defining relative 
clustering strengths for various flow features proved 
difficult, and a damaging lack of robustness is shown 
for three-dimensional leeside flowfields.® 

Mesh adaption on unstructured domains offers a sig- 
nificant benefit over structured mesh adaption — the 
ability to locally insert and delete nodes. Much of 
the research in this area has gone into global remesh- 
ing using isotropic cells B0.BB0.E] Grids composed 
of (nearly) isotropic cells quickly become prohibitively 
large for hypersonic applications, where the capture 
of essentially one-dimensional flow features, such as 
shocks, requires refinement in all three spatial dimen- 
sions. 

More recently, impressive results using anisotropic 
elements have been reported by Habashi et 
a?.®’®’®’®’®® With their approach, all grid 
adaptions are local operations, as opposed to global 
remeshings. Highly-stretched elements are obtained, 
achieved by equating the interpolation error along 
each edge, again using a spring analogy minimiza- 
tion. The clustering is driven by a posteriori error 
estimation based on second derivatives of the solu- 
tion. This sort of clustering is intended to reduce 
the interpolation error in a piecewise-linear data 
representation,®® but is not necessarily driven by 
the flow physics, and can lead to excessive clustering 
or conflicting requirements in certain regions, such as 
a bow shock or stagnation point. Figure |I] presents 
an illustrative pictorial based on the results of Ait- 
Ali-Yahia et aZ,™ where 18 cells were driven into the 
bow shock by gradient based clustering. A shock is 
pictured on the left side of the figure with an ideal 
mesh for a three-point stencil. On the right-hand 
side is a mesh dictated by curvature-based clustering, 
clearly containing more points than necessary. 

Roe®® h as applied the concept of local node move- 
ments to a scalar advection problem using a fluctua- 
tion splitting solver. His analysis reveals that a char- 
acteristic mesh results, with far fewer points required 
than curvature clustering would imply. His method 
is based on the minimization of an objective func- 
tion formed by taking derivatives of the fluctuation 
splitting scheme. A differentiable high-resolution lin- 
ear scheme, which in general cannot be monotonic, 23 
was chosen. For a complex solution field or for sys- 
tems of equations it may not be possible to achieve a 
perfectly-aligned characteristic mesh, in which case the 


non-monotonic property would most likely be detri- 
mental or even fatal to the solution. It is not clear 
how to extend the differentiability requirement to a 
high-resolution, non-linear monotonic scheme. 

The current study seeks to borrow from the full suite 
of aggressive anisotropic unstructured mesh adaptions 
that have been developed for traditional finite vol- 
ume or finite element methods, and to apply them in 
novel ways to the inherently multi-dimensional fluc- 
tuation splitting distribution scheme. Prior work®’® 
has shown fluctuation splitting can deliver greater ac- 
curacy on coarser meshes than finite volume when 
applied to model problems with features representative 
of those present in hypersonic flows. A mesh adap- 
tion strategy developed in conjunction with fluctuation 
splitting that is robust and general yet does not overly 
cluster to shocks could contribute significantly to low- 
ering the cost of computational aerothermodynamics. 

This paper judges the potential of such an adaption 
strategy for a fluctuation splitting solver by compar- 
ing its performance versus a state-of-the-art solution- 
adaptive strategy, based on a posteriori error estimates 
for a finite volume solver, as applied to several scaler 
model problems representative of some of the features 
present in compressible gas dynamics. The govern- 
ing equations are stated along with a brief outline of 
the finite volume and fluctuation splitting discretiza- 
tions. Descriptions of the a posteriori error estimation 
and the basic local adaption techniques are followed 
by details of the adaption strategy developed for the 
present study, constructed as a series of minimization 
operations on the cell fluctuations. Applications to 
linear advection, non-linear advection, and advection- 
diffusion lead to conclusions about the suitability of 
the present method for extension to compressible gas 
dynamics. 

Governing Equations 

The scalar advection-diffusion equation, 

[/ t + VT = V’(/iW) (1) 

is considered in two spatial dimensions on triangu- 
lated unstructured domains. The physical domain is 
chosen to be the unit square in the second quadrant. 
Steady-state solutions are sought through pseudo-time 
relaxation. 

A linear equation is obtained with, 

F = XU (2) 

while non-linear advection is obtained with, 

F={\u\ U), /x = 0 (3) 

Finite Volume 

The finite volume discretization is constructed as an 
edge-based implementation of the approximate Rie- 
mann solver due to Roe. 123 The numerical flux at a 
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median-dual control-volume face is, 

^ [Fi n + F ou ^j -h - \ \A-h\{U out - U in ) (4) 
with the flux Jacobian defined, 


A = 


dF 

dU 


(5) 


Node-based limited reconstruction for second-order ac- 
curacy follows Barth, 23 

U face = U 0 + 4’(VU) 0 -r (6) 


Diffusion terms are discretized as a finite element 
distribution such that the contribution to node i of a 
triangle is given by, 


The absolute values of the advection speeds are limited 
as, 


H if |a| > £( a ) 

yjr if l“l < £ («) 


(14) 


and similarly for |/3|. The small parameters for Eqn. 0 
are obtained as, 


e (a) = x max 


e (/3) = 2 max 


0, l\fi\ ■ (A — Ai), i\n\ ■ (A 2 — A) 
0, ^ 3^3 • (A-2 — A ) , £ 3 fi3 -(A — A 3 ) 


(15) 


The advective fluctuations are distributed to the 
nodes of the triangle as, 


Pi = - ( 7 ) 

3 = 1 

with the convention that edge i is opposite to node i. 
Fluctuation Splitting 

The fluctuation splitting discretization is the multi- 
dimensional upwind residual distribution scheme of 
Sidilkover. 23 The advective fluctuation is defined on a 
triangular cell, 

<!> = -[ y-Fdn = p + p r ' ( 8 ) 

J n 

where, 

P = a(U 1 -U 2 ), r =P(U 3 ~U 2 ) (9) 

and, 

a=^A-n 1 , (3=^A-h 3 (10) 

The fluctuations are limited to achieve linearity preser- 
vation as, 

A = ^ + *(^r) 

f n =4?-i>MS)4>’< (ll) 

and artificial dissipation is introduced for upwind 
monotonicity, 

P* = sign(a)p\ </>'"= sign (/3)^*" (12) 

Expansion shocks are eliminated by extending the 
one-dimensional eigenvalue limiting of Harten and Hy- 
men® to multiple dimensions by searching for expan- 
sions in the £ and directions separately. The artificial 
dissipation can be recast, with Eqn. 0, as, 

p S = \a\(Ui — U 2 ) + sign (a)'4>(j) v 

P* = mU 3 -U 2 )~ sign(/3)^ (13) 


node 1 •*— i (4>* S — p ^ ) 

node 2 <- + p S ) + ^ (</>*" + P") 

node 3 <— — p V ) (16) 

Diffusive terms are computed as a finite element dis- 
cretization and distributed as in Eqn. 0. 

Further details of the specific implementations 
for both the finite volume and fluctuation splitting 
schemes can be found in Ref. G5. 

It will be demonstrated that fluctuation splitting 
can provide exact solutions to multi-dimensional lin- 
ear advection problems, Eqns. [j] and [| with p = 0, 
on a properly aligned mesh. The counter-case, of 
a poorly aligned mesh producing significant artificial 
cross-diffusion, has been covered in Ref. [24|. 

The construction of a properly aligned mesh begins 
by looking at the elemental fluctuation, Eqn. g, which 
can be written, 

<t>=\ [(J7i - UPhfi! + (U 3 - U 2 )i 3 n 3 \-\ (17) 

If one edge of the triangular element, without loss of 
generality say edge 1 (edge 1 is opposite to node 1), 
is parallel to the cellular advection velocity, then 

fii-A = 0, and the fluctuation reduces to, 

(p = - X-h 3 £ 3 (U 3 — U 2 ) (18) 

Now (f> — > 0 as U 2 —■ > U 3 , irrespective of the value 
of U\. Since the artificial dissipation scales propor- 
tionally with the cell fluctuation, the exact isentropic 
solution will be captured as the fluctuation splitting 
scheme converges to U 2 = U 3 . 

A mesh such that each cell has one edge aligned with 
its averaged advection velocity is analogous to a char- 
acteristic mesh. Unfortunately, this edge-alignment 
concept will not, in general, produce an exact solution 
in the presence of physical diffusion because the correct 
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result will typically be U 2 7^ E/3 even if edge 1 is aligned 

with A. Additionally, extensions to systems possessing 
multiple or imaginary characteristics are difficult to 
conceive. 


Curvature Clustering 

Current state of the art for local unstructured 
mesh adaption is based on a posteriori error estima- 
tion 3QEJ 00,0>EllI3II3>f23 The error estimates can 
be derived either by looking at the leading truncation 
error terms in the spatial discretization or by consider- 
ing the solution interpolation error. In practice, either 
approach reduces, for second-order-accurate spatial 
discretizations, to a check on the curvature of the so- 
lution. If isotropic cells are desired, the magnitude of 
the local curvature inversely dictates the element sizes, 
while for anisotropic adaption directional derivatives 
can be used to stretch elements. 

Habashi 1 ^ defines the a posteriori error estimate on 
an edge to vary like, 


\E r \ ~E 2 


d 2 U 
dr 2 


(19) 


or as “the edge length squared times the second deriva- 
tive of the solution.” In the finite volume context, one 
simple and efficient way to construct an edge-based 
error estimate along edge 01 is, 


\E r \ 


VE/r^-VE/o-^ 


toi 


L 01 



(20) 


Barth’s edge-based finite volume scheme, used here, 
already requires the gradient computations, making 
the error estimate a trivial step. 

Having defined an error estimate for all edges, the 
adaption strategy seeks to reduce the magnitude of the 
estimates while anisotropically stretching cells to equi- 
distribute the error across all edges. The adaption is 
performed locally using the four basic operations: edge 
swapping, point deletion, point insertion, and nodal 
displacement. 

One simple method to improve a mesh is to swap 
edges between nodes, altering the local connectivity. 
If the triangles to either side of an edge form a con- 
vex quadrilateral, then that edge is a candidate to be 
swapped. If the error estimate for the swapped connec- 
tivity is smaller than for the current edge connectivity, 
then the edge is swapped. In practice, an error thresh- 
old is employed with all the local operations discussed 
herein to avoid limit cycles in smooth regions of the 
solution. 

If all edges emanating from a node have very small 
error estimates, then that node is deleted. If a given 
edge has an excessively large error, then that edge is 
split by adding a node at the midpoint. 


The approach for moving nodes uses the spring anal- 
ogy, similar to the techniques presented by Gnoffo^ 
and Ait-Ali-Yahia et alF n The springs are taken to 
be the mesh edges. The spring constants are the edge 
error estimates. A minimization is then sought for a 
potential energy formed as, 

E ^ \itfi ( 21 ) 

3 

where j represents all distance-one neighbors to node 0. 
The minimum is obtained by setting the derivative to 
zero, giving the requirement, 

El \Er\j To'j = 0 ( 22 ) 

3 

where O' indicates the new, minimum energy position 
of node 0. Holding the edge error estimate constant 
during the displacement, an updated position vector 
can be obtained, 


E \ Er \i + r 0 'o) = 0 

3 


Zoo' — — H)'o — 


E ? 3 r 0 3 

Yhj \Er\j 


(23) 

(24) 


Equation gives the position vector pointing from 
the current location of node 0 to its adapted loca- 
tion. This adapted location is an attempt to optimally 
equate the scaled error estimates over all edges con- 
nected to the current node. 

Following the recommendations of Dompierre et al 13 
and Habashi et al ^ the nodal displacement technique 
is used as a smoother between applications of point 
insertion/deletion and edge swapping. Coupling be- 
tween the solver and the mesh adaptor is maintained 
by iterating the solver between adaption operations. A 
complete adaption cycle consists of the following steps: 

Solve on initial mesh. 

Swap diagonals and iterate solver. 

Move nodes and iterate solver. 

Insert nodes and iterate solver. 

Move nodes and iterate solver. 

Delete nodes and iterate solver. 

Swap diagonals and iterate solver. 

Move nodes and iterate solver. 

The preceding steps are repeated until convergence of 
the entire process. 


Fluctuation Minimization 

Fluctuation splitting, being a distribution scheme, 
discretizes the partial differential equation into a set 
of algebraic relations for the nodal values of the depen- 
dent variable. Convergence of the solver implies that 
the nodal updates for U are driven to zero. However, 
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the nodal updates are formed as contributions from 
all surrounding cells, and the fluctuations defined on 
these cells are not necessarily reduced toward zero, as 
the possibility exists for positive and negative contri- 
butions from neighboring cells to cancel on summation 
about a node. 

Recall that the fluctuation splitting artificial dissi- 
pation terms, Eqn. 0- scale with the cell’s fluctuation, 
and that for linear advection it is possible to eliminate 
the production of artificial dissipation by using an ap- 
propriately designed mesh. Explicitly adapting a mesh 
to align with characteristic directions proves awkward 
to translate into an algorithm, and offers limited hope 
for extending to diffusive problems or systems. 

An alternative tact is chosen whereby the mesh is 
adapted so that the solver will minimize both the cell 
fluctuations, and hence the artificial dissipation, as 
well as the nodal updates. Since there are on the or- 
der of twice as many cells as nodes, the fluctuations 
are minimized in a least-squares sense. The advantage 
of fluctuation-based, rather than characteristic-based, 
adaption is that the fluctuations remain defined for 
diffusive problems and systems. 

Given a sub-optimal mesh and solution cor- 
rupted with artificial dissipation, the present adaption 
method performs a series of local operations, driven 
by predictions of what the fluctuations will be on the 
modified mesh. The sequence of local optimizations is 
iterated to obtain globally improved solutions. 

For edge swapping, the root mean square (RMS) of 
the fluctuations to either side of the edge is compared 
to the RMS of the fluctuations in the swapped config- 
uration. If the swapped RMS is lower, then the edge 
is swapped. If the fluctuations in all cells surrounding 
a node are very small and the fluctuations will remain 
small in the re-triangulated region without that node, 
then the node is deleted. If the fluctuations in the two 
cells to either side of an edge are very large, then that 
edge is split by adding a node at the midpoint. 

For displacing nodes, a scheme to minimize RMS 
of fluctuations has been presented by Roe.®® The 
development presented by Roe uses the same mini- 
mization scheme to evolve the solution as is used to 
drive the nodal positioning. That type of fluctuation 
splitting solver has a central difference flavor and is 
not monotonic. The present procedure incorporates 
the upwind, non-linear fluctuation splitting algorithm 
of Sidilkover into some of the mesh movement strate- 
gies of Roe. 

At a given node, the nodal displacement is computed 
as a first step and then the solution is updated at the 
new nodal location via a local point-implicit inversion. 
In this manner a global mesh movement sweep can be 
accomplished in conjunction with a single Gaufi-Seidel 
iteration of the solver. In this section the current node 
to be moved is globally numbered node i. Within each 
triangular element the nodes are locally numbered 1-3. 


Derivatives in y are omitted when they exactly follow 
from the derivatives in x. 

The functional to be minimized is defined at the 
node as a sum of contributions from all cells surround- 
ing the node (equivalent to Eqn. 13, p. 248 in Ref. f22j ), 

t * = ^£ s t # ( 25 ) 

T 

for all triangles T containing node i. The weight- 
ing factor, Sy, is a positive scalar, generalizing to a 
symmetric, positive definite matrix for systems. The 
functional is thus a sum of positive semi-definite con- 
tributions from triangles containing the current node. 

The derivative of T with respect to a nodal coordi- 
nate is, 


— = V ( ^ dEj i c a d( t >T \ 
dxi ^ \ 2 dxi T T dxi J 


(26) 


Note that the derivatives in Eqn. ^ represent the 
change in solution values as the discrete mesh is per- 
turbed, and as such are to be interpreted in the context 
of variational calculus, and not as spatial gradients ac- 
cording to the more-familiar multi-variable calculus. 

The minimization of T can be performed using a 
fixed-point iteration to force the derivatives to zero, 


x? +1 


= X, 


dr 

dx^ 



(27) 


Equation can be combined with Eqn. ^ in the form 
of a distribution method of steepest descent, 


A ,: 


<+ dr 


t , dfa_ 

dxi T T dxi 


(28) 


Convergence can be enhanced over the fixed-point 
iteration by using a Newton scheme. Expanding the 
gradient in an approximate Taylor series, 


0 = 


dT 

dxi 


n+1 


dr 

dxi 


A n x, 


d f dr 


An 9 f 9T 

Vl d Vi [dxi 


0 = 


dr 

% 


n+1 


dr 

% 


A" 


dxi V dx 


d (dr 


dxi \dyi 


An d (dr 

Vl dyi \dyt 


(29) 


(30) 


leading to the form, 



o 2 r 

a 2 r 

-1 

f dr ) 

dxf 

dxidyi 

< 

dxi \ 

d 2 r 

a 2 r 


dT I 

dxidyi 

dy 2 


l dy, ) 


(31) 


In Ref. |2T : , Roe suggests neglecting the off-diagonal 
terms in Eqn. m 
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Second derivatives of the objective function are, 


d 2 T 

fa 2 


= £ 


T L 


d 2 Zj 
2 dx 2 


2 <t>T 


_ dZ T dcpT _ f d(j ) T 

."".T 


r^T 


dxi dxi 

d 2 (j)j' 
dx 2 


V 

(32) 


<9 2 T 

dxidyi 



+ S T 


/ ds T 

v&u % 

0£i % 


<95 T d<ft T \ 

dxi dyi ) 

6 ~ -£^L 

T “ T dxi.dy.i_ 


# d 2 ~ T 

2 dxidyi 


(33) 


Explicit expressions for the derivatives of the fluctu- 
ations are grouped in the appendix. For application to 
fluctuation splitting schemes other than Sidilkover’s, 
only those equations in the appendix would change. 

For 5, Roe 22 chooses 5y = for which the deriva- 
tives are, 



dzij 1 dSj 

(34) 


dxi S 2 dxi 

3 2 S t 

2 fdS T V 1 d 2 Sj 

(35) 

dx 2 

^ dxi ) S 2 dx 2 

d 2 Z T 

2 dS T dS T 1 d 2 Sj 

(36) 

dx^yi 

Sj dxi dyt Sj dxidyi 


The area of a triangle is, 

St = ^ [a=i(3/2—2/3) + £2(2/3-221) + £3(2/1 -3/2)] (37) 


leading to the derivatives, 


dS- T _ 1/3 - 2/1 
dx 2 2 


dSj £1 — £3 
dy 2 2 


d 2 S T _ d 2 ^ _ d 2 Sj 
dx 2 dy 2 dx 2 dy 2 


(38) 


(39) 


This choice of weighting emphasizes fluctuations on 
the smaller cells. 

An alternative is to weight all cells equally, with 
St=1. Two other obvious choices for weighting are 
St=St, emphasizing fluctuations on the larger cells, 
and 5 t = Jj-, for even stronger emphasis on the 
smaller cells. The derivatives for this latter case are, 
with Eqn. 


dZ t _ _ 2 dS T 

dxi Sj dxi 

d 2 z T _ 6 /as T \ 2 

dx 2 ~ \ JhTi ) 

d 2 Z T _ 6 dSjdSj 

dx^yi dxi dyi 


(40) 

(41) 

(42) 


A complete adaption cycle follows the same steps as 
outlined for curvature-clustering adaption. 



X 



Fig. 2 Unadapted finite volume results for linear 
advection of shear. Solution contours vary on (0,1) 
with 0.1 increments. 


Results 

Four scalar test cases are considered: uniform ad- 
vection, circular advection, non-linear advection, and 
advection-diffusion. The emphasis of the current study 
is to evaluate the performance of the adaption schemes 
using only a few adaption cycles, typically 1-4, on 
each case. In this manner, significant improvements in 
the solutions are sought with little additional overhead 
from the adaptions. It should be noted that continued 
application of either adaption strategy would lead to 
continually refined solutions, though at the cost of ever 
increasing numbers of nodes. 

Linear Advection 

The first test case is for linear advection of a shear 
discontinuity at 45°, A=(l,l). Inflow conditions are 
U(—l,y) = 0, U(x, 0) = 1. The starting 121-node 
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Fig. 3 Adapted mesh with greatly improved finite 
volume solution for 45 ° shear. Solution contours 
vary on (0,1) with 0.1 increments. 

mesh and converged finite volume solution, using the 
compressive van Albada limiter, are shown in Fig. 
Excessive spreading of the contour lines, indicative of 
the magnitude of non-physical dissipation corrupting 
the solution, is seen. The fluctuation splitting solution 
on this mesh, using the Minmod limiter, is similar in 
character, though somewhat less diffusive, 1111 to the fi- 
nite volume solution. 

Four cycles of the curvature-clustering adaption 
with the finite volume solver result in the 103-node 
mesh and solution of Fig. 0. While the adapted mesh is 
highly anisotropic, a significant improvement in reso- 
lution of the shear has been obtained with a 15 percent 
reduction in the number of nodes. 

However, a single application of only the edge- 
swapping and point-deletion routines with the fluc- 
tuation splitting solver results in the optimal mesh 
and exact solution, Fig. f|, obtained with no interior 
nodes. The six retained boundary nodes are a 95 per- 
cent reduction in the number of nodes. Applying the 
finite volume solver to this optimal fluctuation split- 
ting mesh gives the diffused result of Fig. 0, showing 



1 T 

0.8 - /W 

Y 

0.6 - /w' 

0.4 - 

0.2 - Aw 

0 L l l l l l 

-1 - 0.8 - 0.6 - 0.4 - 0.2 0 

X 

Fig. 4 Optimal fluctuation splitting adaption and 
exact solution for 45 ° shear case. 



ol 1 1 1 1 1 

-1 - 0.8 - 0.6 - 0.4 - 0.2 0 


X 

Fig. 5 Finite volume solution on the optimal fluc- 
tuation splitting mesh. 
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Fig. 6 Starting mesh and converged finite volume 
solution, A = (y, —x). Solution contours vary on 
(0,1) with 0.1 increments. 

very poor containment of the leading and trailing con- 
tour levels. 

Circular Advection 

The second case considered is for circular advection, 
with a variable advection velocity of A = (y, — x). The 
inflow profile is U(x, 0) = 0 on x = [—1,0.7] and 
t/(x, 0) = 1 on x — [—0.6,0]. The initial mesh and 
highly diffused finite volume solution on this mesh are 
shown in Fig. 0. As before, the fluctuation splitting so- 
lution on this mesh is similar in character, though less 
diffusive than the finite volume result. Three adap- 
tion cycles of the curvature clustering with the finite 
volume solver are applied, leading to the mesh and so- 
lution of Fig. 0. The number of nodes has increased 
from 121 to 146, but the solution has been dramati- 
cally improved. 

For this problem, an optimal mesh can be con- 
structed following the characteristic-alignment guide- 
lines. Such a mesh is shown in Fig. 0, along with the 
exact solution as computed by the fluctuation split- 
ting solver. This mesh contains only 10 nodes. The 


Fig. 7 Converged finite volume solution and 
curvature-clustered mesh after three grid-adaption 
cycles, A = (y,—x). Solution contours vary on (0,1) 
with 0.1 increments. 

finite volume solution on this mesh, Fig. 0. is clearly 
far inferior to the exact fluctuation splitting solution. 

In practice, though, the local adaption algorithm 
for fluctuation splitting presented here is not able to 
reduce the mesh to the ‘optimal’ grid for this case. 
This is because each local operation is constrained 
to be in the direction of a local improvement only. 
In order to reduce the starting grid to the ‘optimal’ 
mesh, a global communication is required between 
nodes to know that a local move in the non-optimal 
direction will be counter-acted by changes occurring at 
another node. The present method is applied to the 
circular advection problem for one cycle, using four 
sub-iterates of the diagonalized Newton scheme dur- 
ing the nodal displacement step. The resulting mesh, 
containing 70 nodes, and fluctuation splitting solution 
are shown in Fig. |TIJ The fluctuation splitting result 
is of comparable accuracy, but rougher at the edges 
of the shear, to the adapted finite volume solution ob- 
tained on a grid containing twice as many points. 


8 of EZ3 


American Institute of Aeronautics and Astronautics Paper 99-3254 




Fig. 8 Optimized mesh, created by hand, and ex- 
act solution to circular advection problem, using 
fluctuation splitting. 



Fig. 9 Finite volume solution on the optimal fluc- 
tuation splitting mesh for circular advection. 




Fig. 10 Mesh and fluctuation splitting solution af- 
ter one adaption cycle for circular advection. Con- 
tours vary on (0,1), with 0.1 increment. 

Non-linear Advection 

A non-linear advection case is constructed contain- 
ing a symmetric compression fan that coalesces into a 
vertical shock at (— |, I). Centered expansion fans sit 
at (-1,0) and (0,0). The inflow profile is, 

U(x,0) = — 2x — 1 on a; = (—1,0) 

f7(— 1,0) = 17(0,0) = 0 (43) 

The unadapted mesh and extremely diffused finite vol- 
ume solution are shown in Fig. [TT]. Three full cycles 
of the curvature-clustering adaption nearly double the 
number of nodes to 237, yielding the mesh and solution 
in Fig. 0 The solution on this mesh shows a dramatic 
improvement for shock thickness, shock speed, point of 
coalescence, and preservation of extremum in smooth 
regions. Note, however, that there is some asymmetry 
between the expansion fans toward the top of the do- 
main, that the shock is not entirely straight, and that 
the compression fan begins to coalesce into a shock at 
y = 0.4, instead of at y = 0.5, the correct location. 

The fluctuation splitting adaption scheme is also 
applied for three cycles to the same problem and ini- 


9 of DU 


American Institute of Aeronautics and Astronautics Paper 99-3254 



Fig. 11 Initial mesh and finite volume solution for 
non-linear advection case. Contours vary on (-1,1), 
with 0.1 increment. 

tial grid. The adapted mesh, containing 206 nodes, 
and corresponding solution are shown in Fig. 0. For 
this test case, there are pros and cons for both sets 
of results. The adapted fluctuation splitting solution, 
using 14 percent fewer nodes, exhibits slightly greater 
accuracy than the adapted finite volume solution, par- 
ticularly in the expansion fan symmetry and shock 
coalescence point. One feature that is better resolved 
in the finite volume solution is the extremum between 
the compression and expansion fans on the lower right- 
hand side. The fluctuation splitting shock is broader 
near the coalescence point but tapers to a comparable 
crispness with the finite volume result, and is a little 
straighter. The fluctuation splitting shock at the out- 
flow is ever so slightly offset to the right from x = —0.5, 
the correct location. 

The finite volume and fluctuation splitting meshes, 
Figs. 0 and 0 have similar character in the expan- 
sion and compression fans. The fluctuation splitting 
adaption does not cluster as many points to the shock, 
resulting in overall 14 percent fewer nodes. This fact 
that the fluctuation splitting adaption scheme can re- 


Fig. 12 Curvature-clustered mesh and finite vol- 
ume solution after three adaption cycles, non-linear 
advection case. Contours vary on (-1,1), with 0.1 
increment. 

solve crisp shocks without excessive clustering normal 
to the shock may have favorable implications when 
considering fluid dynamic problems with extremely 
strong shocks in the vicinity of more subtle, though 
still important, features, such as an entropy layer. 

Advection-Diffusion 

The final case is for an advection-diffusion problem 
due to Smith and Hutton. 52 The advection velocity is, 

A= (2y(l-x 2 ), -2x(l-y 2 )) (44) 

The inflow profile is, 

U(x, 0) = 1 + tanh(20x + 10) (45) 

The diffusion coefficient is chosen to be a small con- 
stant, fi = 10 — 3 , to be representative of a high- 
Reynolds-number shear. The starting mesh for this 
case is an unstructured isotropic mesh containing 1928 
nodes. The reference solution for this case is taken 
from a fully grid-converged solution on an isotropic 
mesh (20,000 nodes). Both the finite volume and 
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Fig. 13 Mesh and fluctuation splitting solution 
after three adaption cycles, non-linear advection 
case. Contours vary on (-1,1), with 0.1 increment. 

fluctuation splitting solvers are run on the unadapted 
mesh, and the RMS difference in the outflow profiles 
between these solutions on the starting mesh and the 
fully-resolved solution are compared in table [I]. As 
usual, the fluctuation splitting solution is seen to be 
more accurate than the finite volume solver on the un- 
adapted mesh. 

Two cycles of curvature-clustering adaption are ap- 
plied, reducing the number of nodes by two-thirds to 
619, while improving the RMS outflow resolution by 
28 percent, also listed in table |T|. One cycle of the 
fluctuation splitting adaption reduces the number of 
nodes to 695, while still producing some (7 percent) 
improvement in accuracy. 

Concluding Remarks 

Current state of the art for local anisotropic un- 
structured mesh adaption based on a posteriori er- 
ror estimates has been implemented in an edge-based 
structure in conjunction with a finite volume solver. 
This type of adaption results in meshes where the node 


Table 1 RMS difference of solution outflow profile 
relative to reference solution. 



Finite volume 

Fluctuation splitting 

unadapted 

0.0288 

0.0068 

adapted 

0.0208 

0.0063 


densities are clustered to regions of high curvature in 
the solution. Significant improvement in solution ac- 
curacy is verified using this technique on scalar model 
problems. 

Recognizing the remarkable property of the dis- 
cretized fluctuation splitting scheme that multi- 
dimensional advection can be solved exactly when 
one edge of each cell is aligned with the character- 
istic direction, a different mesh adaption scheme is 
proposed. While retaining the mechanics of perform- 
ing only local operations, i. e. point insertion/deletion, 
edge swapping, and nodal displacement, a solution- 
predictive approach is chosen in favor of a posteriori 
curvature clustering. The concept of aligning cell 
edges with characteristic directions is generalized as 
a minimization procedure to allow extension to diffu- 
sion problems and systems. Extending this process 
to non-linear problems led to an implementation of 
eigenvalue limiting for multi-dimensional upwind dis- 
tribution schemes. 

It is seen that, while performing a series of local 
optimizations does lead to globally improved solution 
accuracy and reduced grid sizes, in general a truly 
globally ‘optimal’ mesh is not achieved in a small 
number of adaption cycles. However, the solution- 
predictive adaption in conjunction with the fluctuation 
splitting scheme does provide moderately more accu- 
rate solutions on smaller meshes for comparable num- 
ber of adaption cycles versus the finite volume solver 
with adaption driven by error estimates. 

Considering extensions to three-dimensional hyper- 
sonic flow applications, perhaps the most promising 
difference between the two adaption strategies lies in 
their treatment of shocks. For the a posteriori adap- 
tion, the number of nodes clustered to the shock grows 
as the shock strength grows, which can lead to a bow 
shock dominating the adaption for hypersonic prob- 
lems. In contrast, the minimization of fluctuations 
tends to merely align the grid with the shock, leav- 
ing the points outside the shock largely unaffected. 

Appendix 

The variations of the cell fluctuations and the de- 
pendent variable nodal values with respect to changes 
in the nodal location are developed for the advec- 
tive distribution scheme of Sidilkover and the diffusive 
Galerkin distribution. 
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Advective Fluctuations 

The advective fluctuation can be manipulated from 
Eqn. |TT| as, 

<t>= \ [(^2 - U 3 )( Vl , -xi) + ( U 3 - Ui)(y 2l - x 2 ) 

+ {U 1 -U 2 ){v 3 ,-x 3 )]-X (46) 
The spatial derivatives take the form, 


d(j) 1 
dx 2 2 
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3 = 1 
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3 = 1 
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d 2 (j) 

dx 2 dy 2 
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(53) 
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(48) 


The derivatives of the cell-averaged flux Jacobian 
will depend upon the particular flux function. How- 
ever, since A is a weighted average of three nodal 
values, 


dA 

dx 2 


1 dA 2 
3 dx 2 


(49) 


Second derivatives of A are developed in an analogous 
manner to its first derivatives. 

Diffusive Fluctuations 

Derivatives of the diffusive fluctuations take the 
form, 


dx 2 


4SV 


/1 dfi 

1 £>s t \ 

\ftdx 2 

St d x 2 / 


3 = 1 


,dU- 


+ ^2 -Q^ 1" (®3 — Xi)(Ui — U 3 ) 


(54) 


For uniform advection, = 0. For circular advec- 


tion, = (0,-4) and = (§,0). For non-linear 
problems, in general, 


dA 

dx 2 


1 dU 2 
3 dx 2 


(50) 


In keeping with the GauB-Seidel update philosophy, 
only the derivative of the solution value at the current 
node is retained in the last term of Eqns. S3 and [TH| . 
That is, if the current node is designated node 2 of 
the triangle under consideration, then is 

retained while ~ ~ 0 is assumed. 

Second derivatives of the advective fluctuation with 
respect to variation of a nodal location follow from 
Eqns. S3 and g8|, incorporating the approximation 

dUi ^ dUz ^ q 


dx2 


dx2 


d 2 (f _ tTT TT ,dA* , „ „ an du 2 , h „ fd 2 u 2 

~k ~ 2 — {U 1 — U 3 )— hr 2 n 2 -^ — w 1- —n 2 -A 2 

<7X2 <7X2 OX 2 2 ( 7 X 2 

+ (51) 


j'= 1 
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9y 2 


45 t 


/ 1 S/r 

1 ds T \ 

\ydy 2 

St <9y 2 / 


i=i 


+ ^^ + (ys-yi)(Ui-u 3 ) 

dy 2 


(55) 


where the assumption ~ 0 has already 

been applied. The second derivatives of the diffusive 
distribution follow as, 


92 

dxl 


fif 2 d 2 U 2 ( y dSj 1 dfi 


4St dx\ 


\4Sf dx 2 457 dx 2 


dU 

tl-^ + {x 3 -xi){ u i- u 3 ) 


dfi N 2 


4/tSt \dx 2 
2 " 


1 d 2 y fj, d 2 Sj fj, f <9St\ 
” 4SV dx% + 4 S 2 dxl ~ ) 
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g W2 n j .n 2+ ^-— j— (56) 
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d 2 p 2 y £ 2 d 2 U2 ( y dSj 1 dy \ 
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(58) 


is not used in the calculation of the solution, but only 
employed to provide an estimate of the solution vari- 
ation with respect to nodal displacements, providing 
the forcing functions for the mesh adaption. 

Assembling the linear upwind advective distribu- 
tions with the diffusive contributions from the sur- 
rounding cells and solving for the steady state value 
of the current node yields, with U 2 = Ui, 


UiY,( a+ -F + ^=Y,[ a+ Ui-p-Us 

- (Uihni + U 3 l 3 n 3 )-rii 


(60) 


where, 


. adzlal 1 d= sign(a) 

a = — = a 

2 2 

= (61) 

The variation of the nodal solution with respect to 
nodal displacement can now be evaluated as, 


The derivatives of the cell-averaged diffusion coeffi- 
cient scale like, 

d 2 / 2 1 d 2 y 2 , 

dx 2 3 dx 2 dx 2 3 dx 2 

Dependent Variable 

Evaluating directly from the high-resolution 
non-linear fluctuation splitting scheme is impractical. 
Limiters such as Minmod are not continuously differ- 
entiable, and while the van Albada limiter is differen- 
tiable its use does not lead to a convenient explicit 
form from which to evaluate As an approxi- 

mation to (jf for the non-linear scheme, derivatives 
are sought using linear distribution schemes. Two 
linear choices for fluctuation splitting are to use a 
linearity preserving (second-order spatial accuracy), 
non-monotonic distribution or an upwind, monotonic 
first-order distribution. The linearity-preserving, non- 
monotonic scheme is of the Lax-Wendroff type,® 1 and 
tends to produce dispersion waves in response to nodal 
displacements. This behavior tends to under-predict 
the change in solution value at the perturbed node 
relative to the fluctuation splitting scheme employed 
herein. 

The linear upwind scheme is obtained from the 
present fluctuation splitting scheme by discarding the 
limiter. This scheme exhibits a dependency on the 
numbering of nodes within a triangle. To alleviate 
this dependency, the approximation to is built by 
looping over all cells connected to node i and locally 
renumbering the nodes within each triangle so that 
the current node is designated as node 2 of the trian- 
gle. It is emphasized that the linear upwind scheme 
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Recall that the solution at surrounding nodes is held 
fixed during displacements of the current node. 

The derivatives of the functions in Eqn. m are de- 
fined, 
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dx 1 if, a>0 dx \ 0 , /3 > 0 


Further, Eqn. E3 leads to, 


da 

& 

1„ „ 

dA 

da 

A x 

1„ „ 

dA 

(64) 

dx 2 

“T 

+ 2^i n i ' 

dx 2 ’ 

dy 2 

2~ 

+ 2^ in i ' 

dy 2 

dp 


1„ 

dA 

dp 

A x 

1„ , 

dA 

(65) 

dx 2 

~~2 

~2 W 

dx 2 ’ 

dy2 

~T 

~2 W 

dy2 


For non-linear problems a circular definition arises 
for jSy and One option is to neglect the deriva- 
tives of A in Eqns. and ^s|. Another option would 
be to lag these same terms from the previous iteration 
level. 

The second derivatives of the solution at the current 
node with respect to variations in the position of the 
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current node can be developed from the first derivative 

expression in Eqn. again with the approximation 
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The remaining derivatives to be specified follow from 
Eqns. [TT] and [k»|. 
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